#include "StdAfx.h"

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include "ephcom.h"
#include "EphemDE.h"

#define _CRT_SECURE_NO_WARNINGS
#pragma warning(disable : 4996)

#ifndef _ASSERT
#define  _ASSERT
#endif

/*
	Get the next group header in the file.
	Group header lines are 12 characters long, of the form "GROUP   nnnn".

	Parameters:
	group - the GROUP header we read
	expected - the header we expected
	infile - the file pointer to read
*/
char *ephcom_nxtgrp(char *group, char *expected, FILE *infile) 
{
	char readbuf[EPHCOM_MAXLINE + 1];

	fgets(readbuf, EPHCOM_MAXLINE, infile); /* Blank Line     */
	fgets(readbuf, EPHCOM_MAXLINE, infile); /* "GROUP   dddd\n" */
	strncpy(group, readbuf, 12);
	group[12] = '\0';
	if (strncmp(group, expected, 12) != 0) 
	{
		fprintf(stderr, "Badly formed header; \"%s\" not found.\n\n", expected);
		exit(1);
	}
	fgets(readbuf, EPHCOM_MAXLINE, infile); /* Blank Line    */

	return(NULL);
}


/*
   Read a JPL Ephemeris ASCII header from the file pointed to by infp
   and store values in header structure.  Write any errors to stderr.
*/
int ephcom_readascii_header(FILE *infp, struct ephcom_Header *header) 
{
	char group[13]; /* To store the "GROUP" header line information */
	int i, j;
	int iword; /* word number we're reading in a line */
	int blockbytes; /* number of bytes in a block, equals 8 * ncoeff */

	char readbuf[EPHCOM_MAXLINE + 1];

	/*
	First header line: KSIZE= # NCOEFF= #
	*/
	if (infp != stdin) 
		rewind(infp);
	fgets(readbuf, EPHCOM_MAXLINE, infp);
	sscanf(readbuf, "%*6s%6d%*11s%6d", &header->ksize, &header->ncoeff);
	blockbytes = 8 * header->ncoeff;  /* The size of a double, times # of doubles/block */
	if (header->ksize != 2*header->ncoeff) 
	{
		fprintf(stderr, "Badly formed header; KSIZE <> 2*NCOEFF\n\n");
		exit(1);
	}
	/*
	GROUP 1010: Title of ephemeris (DE/LE number, start JD, end JD)
	*/
	ephcom_nxtgrp(group, "GROUP   1010", infp);
	fgets(header->ttl[0], EPHCOM_MAXLINE, infp);  /* JPL Ephemeris title line */
	if (strncmp(header->ttl[0], "JPL ", 4) != 0) {
		fprintf(stderr,"\nERROR: file is not a JPL ASCII header file.\n\n");
		exit(1);
	}
	fgets(header->ttl[1], EPHCOM_MAXLINE, infp);  /* Start epoch */
	fgets(header->ttl[2], EPHCOM_MAXLINE, infp);  /* Finish epoch */
	/*
	Convert any newlines or tabs to single spaces.
	*/
	for (i=0; i<3; i++) 
	{
		for (j=0; j<84; j++)
			if (isspace(header->ttl[i][j]))  header->ttl[i][j] = ' ';
		header->ttl[i][84] = '\0';
	}
	/*
	GROUP 1030: Start and End JD, timestep (in JD) per block.
	*/
	ephcom_nxtgrp(group, "GROUP   1030", infp);
	fgets(readbuf, EPHCOM_MAXLINE, infp);
	sscanf(readbuf, " %lE %lE %lE", &header->ss[0], &header->ss[1], &header->ss[2]);
	/*
	GROUP 1040: Constant names.
	*/
	ephcom_nxtgrp(group, "GROUP   1040", infp);
	fgets(readbuf, EPHCOM_MAXLINE, infp);
	header->ncon = atoi(readbuf);
	/*
	Now read the constant names, 10 per line, each 6 characters long
	preceded by 2 blanks.  Pad names with blanks to make 6 characters.
	*/
	for (i=0; i<header->ncon; ) 
	{
		fgets(readbuf, EPHCOM_MAXLINE, infp);
		if ((j = strlen(readbuf)) < 81) 
		{  /* Pad end with blanks for copying */
			while (j < 81) readbuf[j] = ' ';
		}
		for (iword=0; iword<10 && i<header->ncon; iword++, i++) 
		{
			strncpy(header->cnam[i], &readbuf[2 + iword*8], 6);
			header->cnam[i][6] = '\0';
		}
	}
	/*
	GROUP 1041: Constant values.
	*/
	ephcom_nxtgrp(group, "GROUP   1041", infp);
	fgets(readbuf, EPHCOM_MAXLINE, infp);
	header->nval = atoi(readbuf);
	if (header->nval != header->ncon)
	{
		fprintf(stderr,"Error: number of constants and values not equal.\n\n");
		exit(1);
	}
	/*
	Now read constant values, 3 per line, 26 characters each.
	*/
	for (i = 0; i < header->ncon; i += 3) 
	{
		fgets(readbuf, EPHCOM_MAXLINE, infp);
		for (j=0; j< (int)strlen(readbuf); j++)
			if (tolower(readbuf[j]) == 'd') 
				readbuf[j] = 'E'; /* exponent is 'E' */
		sscanf(readbuf, "%lE %lE %lE",
			&header->cval[i], &header->cval[i+1], &header->cval[i+2]);
	}
	/*
	GROUP 1050: Constant values.
	*/
	ephcom_nxtgrp(group, "GROUP   1050", infp);
	for (i =0; i < 3; i++) 
	{
		fgets(readbuf, EPHCOM_MAXLINE, infp); /* Read line of 13 6-digit integers */
		for (j = 0; j < 12; j++) 
		{
			header->ipt[j][i] = atoi(&readbuf[6*j]);
		}
		header->lpt[i] = atoi(&readbuf[6*12]);
	}
	/*
	If there are no coefficients for an ipt[i][] object (i.e., ipt[i][1]==0),
	then ipt[i][0] should contain the value of the next available coefficient
	number rather than 0, as per communication of Myles Standish to Paul Hardy
	on preferred format of ephemeris headers.

	If there are no libration coefficients (i.e., lpt[1]==0), then lpt[0]
	should contain the value of the next available coefficient number rather
	than 0 as well, as per the same communication from Myles Standish.
	*/
	/* First set j to maximum index into ipt[] that has coefficients */
	j = 0;
	for (i=1; i<12; i++)
	{
		if (header->ipt[i][1] > 0 && header->ipt[i][0] > j)
			j = i;
	}
	/* Now set j to next available index count. */
	if (header->lpt[1] > 0 && header->lpt[0] > j)
		j = header->lpt[1] + header->lpt[1] * header->lpt[2] * 3;
	else
		j = header->ipt[j][0] +	header->ipt[j][1] * header->ipt[j][2] * (j==11 ? 2 : 3);

	for (i=1; i<12; i++)
	{
		if (header->ipt[i][0] == 0) 
			header->ipt[i][0] = j;
	}
	if (header->lpt[0] == 0) 
		header->lpt[0] = j;
	/*
	Set the maximum number of Chebyshev coefficients possible for this file,
	to initialize position and velocity Chebyshev coefficient arrays during
	Chebyshev interpolation.
	*/
	header->maxcheby = 0;
	for (i=0; i<12; i++)
		if (header->ipt[i][1] > header->maxcheby)
			header->maxcheby = header->ipt[i][1];
	if (header->lpt[1] > header->maxcheby)
		header->maxcheby = header->lpt[1];

	header->au = 0.0;
	header->emrat = 0.0;
	header->numde = 0;
	for (i = 0; i < header->ncon; i++) 
	{
		if (strncmp(header->cnam[i], "AU    ", 6) == 0)
			header->au = header->cval[i];
		else if (strncmp(header->cnam[i], "EMRAT ", 6) == 0)
			header->emrat = header->cval[i];
		else if (strncmp(header->cnam[i], "DENUM ", 6) == 0)
			header->numde = (int)header->cval[i];
		else if (strncmp(header->cnam[i], "CLIGHT", 6) == 0)
			header->clight = header->cval[i];
		else if (strncmp(header->cnam[i], "LENUM ", 6) == 0)
			header->numle = (int)header->cval[i];
	}
	if (header->numle == 0) header->numle = header->numde;
	/*
	GROUP 1070: Constant values.
	*/
	ephcom_nxtgrp(group, "GROUP   1070", infp);
	/*
	Now we're pointing to the first block of coefficient data, after header.
	Return at the point where we can start reading coefficients.
	*/
	return(0);
}

/*
Read a block of data coefficients from a JPL ASCII ephemeris file.
Returns number of coefficients read, 0 at EOF.
*/
int ephcom_readascii_block( FILE *infp, struct ephcom_Header *header, double *datablock) 
{
	int i, j;
	int datalines; /* lines of data we've read */
	int datapoints; /* points of data we've read/converted/written */
	char readbuf[EPHCOM_MAXLINE + 1];
	double val1, val2, val3; /* To read text line with 3 double precision words */

	/*
	First line in an ASCII block will be the block number, followed by
	the number of coefficients.
	*/
	datalines = 0;  /* Not reported, but leave in for debugging */
	datapoints = 0;
	if (fgets(readbuf, EPHCOM_MAXLINE, infp) && !feof(infp)) 
	{
		sscanf(readbuf, "%d %d", &i, &j);
		if (j != header->ncoeff) 
		{
			fprintf(stderr,
				"\nERROR: ASCII data file's %d coefficients/block\n", j);
			fprintf(stderr,
				"       doesn't match ASCII header's %d coefficients/block.\n\n",
				header->ncoeff);
			exit(1);
		}
		datalines++;
		for (i=0; i < header->ncoeff && !feof(infp); i += 3)
		{
			fgets(readbuf, EPHCOM_MAXLINE, infp);
			for (j=0; j< (int)strlen(readbuf); j++)
				if (tolower(readbuf[j]) == 'd')  readbuf[j] = 'e';
			datalines++;
			/*
			This is horrible, but use "%le" here and "%lE in the other
			ASCII data routine (ephcom_readascii_header) so gcc won't try
			to store the formats in the same location and write to them.
			(Problem with gcc not acting like K&R without -traditional flag
			and without -fwritable-strings flag.)
			*/
			sscanf(readbuf, " %le %le %le", &val1, &val2, &val3);
			datablock[i] = val1;
			datapoints++;
			if ((i+1) < header->ncoeff) 
			{
				datablock[i+1] = val2;
				datapoints++;
				if ((i+2) < header->ncoeff) 
				{
					datablock[i+2] = val3;
					datapoints++;
				}
			}
		}
	}
	return(datapoints);
}



/*
Read a JPL Ephemeris data block in binary format.

This is the only routine in this library that accesses a file
as a direct access file, with a specified block number.  The
block number ranges from 0 on up (starting at first data block,
after the 2 header blocks).  Returns the number of coefficients
read, or 0 at EOF.
*/
int ephcom_readbinary_block(
		FILE *infp,                    /* File pointer for direct access file */
		struct ephcom_Header *header,  /* header struct, already filled in    */
		int blocknum,                  /* Data block number, starting with 0  */
		double *datablock              /* returned coefficient data block     */) 
{
	int i;
	long filebyte;

	filebyte = (blocknum + 2) * header->ncoeff * 8; /* 8 bytes per coefficient */
	fseek(infp, filebyte, SEEK_SET);
	// printf("Blocknum: %d, Byte: %d\n", blocknum, filebyte);
#if 0
	for (i=0; !feof(infp) && i<header->ncoeff; i++) 
	{
		datablock[i] = ephcom_indouble(infp);
	}
#else
	i = fread(datablock, sizeof(datablock[0]), header->ncoeff, infp);
#endif
	if (i < header->ncoeff && feof(infp)) 
		i = -1; /* 0 --> EOF */
	return(i); /* Number of coefficients successfuly read (all or nohing). */
}

/*
ephcom_jd2cal() - convert Julian Day to calendar date and time.

tjd: double precision Julian Day
idate: integer year, month, day, hour, minute, second of tjd
calendar_type: -1=Julian; 0=Automatic; 1=Gregorian

If automatic, use Julian calendar for dates before 15 October 1582.

From pp. 604, 606 in the Explanatory Supplement to the Astronomical Almanac.
*/
int ephcom_jd2cal(double tjd, int idate[6], int calendar_type) 
{
	int ihour, imin, isec;
	int j;
	/* From Explanatory Supplement to Astronomical Almanac, pp. 604, 606 */
	int I, J, K, L, N, D, M, Y;

	tjd += 0.5 + 0.5 / 86400.0; /* Round to nearest second */
	j = (int)tjd;  /* Integer Julian Day */
	tjd = (tjd - j) * 24.0;
	ihour = (int)tjd;
	tjd = (tjd - ihour) * 60.0;
	imin = (int)tjd;
	tjd = (tjd - imin) * 60.0;
	isec = (int)tjd;
	/*
	Julian calendar.  Explanatory Supplement to Astronomical Alamanac, p. 606.
	If automatic, use Julian calendar for dates before 15 October 1582.
	*/
	if (calendar_type == -1 || (calendar_type == 0 && j <= 2299160)) 
	{
		J = j + 1402;
		K = (J - 1) / 1461;
		L = J - 1461 * K;
		N = (L - 1) / 365 - L / 1461;
		I = L - 365 * N + 30;
		J = (80 * I) / 2447;
		D = I - (2447 * J) / 80;
		I = J / 11;
		M = J + 2 - 12 * I;
		Y = 4 * K + N + I - 4716;
	}
	/*
	Gregorian calendar.
	*/
	else  /* Explanatory Supplement to Astronomical Almanac, p. 604 */
	{	
		L = j + 68569;
		N = (4 * L) / 146097;
		L = L - (146097 * N + 3) / 4;
		I = (4000 * (L + 1)) / 1461001;
		L = L - (1461 * I) / 4 + 31;
		J = (80 * L) / 2447;
		D = L - (2447 * J) / 80;
		L = J / 11;
		M = J + 2 - 12 * L;
		Y = 100 * (N - 49) + I + L;
	}

	idate[0] = Y;
	idate[1] = M;
	idate[2] = D;
	idate[3] = ihour;
	idate[4] = imin;
	idate[5] = isec;

	return(0);
}


/*
	ephcom_doublstrc2f() - function to convert a string with a double precision
	value written in C to a double precision value that
	FORTRAN creates.  Conversion happens in place.
*/

int ephcom_doublestrc2f(char *buf) 
{
	int istart, istop, exp, edigits;
	double x;

	for (istart=0; isspace(buf[istart]); istart++);
	x = atof(&buf[istart]);
	// printf("x = %E\n", x);
	for (istop=istart; toupper(buf[istop]) != 'E'; istop++);
	exp = atoi(&buf[istop+1]);
	exp++;
	if (exp < 0) 
	{
		buf[istop+2] = '-';
		exp = -exp;
	}
	else 
	{
		buf[istop+2] = '+';
	}
	if (x == 0.0) exp=0;
	if (exp < 100) edigits = 2;
	else if (exp < 1000) edigits = 3;
	else edigits = 4;
	// printf("istart=%d, istop=%d, exp=%d, edigits=%d\n",
	//        istart, istop, exp, edigits);

	while (edigits > 0) {
		buf[istop + edigits + 2] = exp % 10 + '0';
		exp /= 10;
		edigits--;
	}

	buf[istop+1] = 'D';

	while (istop > istart && buf[istop-1] != '.') {
		buf[istop] = buf[istop-1];
		istop--;
	}

	buf[istop] = buf[istop-2];  /* buf[istop-1] == '.' */
	buf[istop-2] = '0';         /* leading zero */

	return(0);
}

/*
	Write header information in ASCII format.
*/
int ephcom_writeascii_header(FILE *outfp, struct ephcom_Header *header) 
{
	double val1, val2, val3; /* To read text line with 3 double precision words */
	int i, j, k;
	int blockbytes; /* number of bytes in a block, equals 8 * ncoeff */
	char spaces[84]="                                                                                \r\n";
	int idate[6];
	char *month[12] = {"JAN","FEB","MAR","APR","MAY","JUN",
		"JUL","AUG","SEP","OCT","NOV","DEC"};

	char writebuf[EPHCOM_MAXLINE + 1];


	/*
	First header line: KSIZE= # NCOEFF= #
	*/
	blockbytes = 8 * header->ncoeff;  /* sizeof(double) * # of doubles/block */
	sprintf(writebuf, "KSIZE=%6d    NCOEFF=%6d", header->ksize, header->ncoeff);
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]);
	fprintf(outfp, writebuf);
	if (header->ksize != 2*header->ncoeff) {
		fprintf(stderr, "Badly formed header; KSIZE <> 2*NCOEFF\r\n\r\n");
		exit(1);
	}
	/*
	GROUP 1010: Title of ephemeris (DE/LE number, start JD, end JD)
	*/
	fprintf(outfp, spaces); /* blank line */
	sprintf(writebuf, "GROUP   1010");
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]);
	fprintf(outfp, writebuf);
	fprintf(outfp, spaces); /* blank line */
	/*
	Header title lines with dates, for example:

	JPL Planetary Ephemeris DE405/LE405
	Start Epoch: JED=  2305424.5 1599 DEC 09 00:00:00
	Final Epoch: JED=  2525008.5 2201 FEB 20 00:00:00
	*/
	sprintf(header->ttl[0],"JPL Planetary Ephemeris DE%03d/LE%03d",
		header->numde, header->numle);
	k = strlen(header->ttl[0]);
	strcpy(&header->ttl[0][k], &spaces[k]);
	ephcom_jd2cal(header->ss[0], idate, 0);
	sprintf(header->ttl[1],"Start Epoch: JED=%11.1f%5d %3s %02d %02d:%02d:%02d",
		header->ss[0], idate[0], month[idate[1]-1], idate[2],
		idate[3], idate[4], idate[5]);
	k = strlen(header->ttl[1]);
	strcpy(&header->ttl[1][k], &spaces[k]);
	ephcom_jd2cal(header->ss[1], idate, 0);
	sprintf(header->ttl[2],"Final Epoch: JED=%11.1f%5d %3s %02d %02d:%02d:%02d",
		header->ss[1], idate[0], month[idate[1]-1], idate[2],
		idate[3], idate[4], idate[5]);
	k = strlen(header->ttl[2]);
	strcpy(&header->ttl[2][k], &spaces[k]);

	/*
	Don't print trailing blanks at the end of these 3 lines.
	*/
	for (i=0; i<3; i++) {
		// strncpy(writebuf, header->ttl[i], 85);
		// for (j=83; isspace(writebuf[j]) || writebuf[j]=='\0'; j--) writebuf[j]='\0';
		// if (i > 0) writebuf[++j] = ' '; /* To match end space in JPL header */
		strncpy(writebuf, header->ttl[i], 80);
		writebuf[80] = '\r';
		writebuf[81] = '\n';
		writebuf[82] = '\0';
		fprintf(outfp, writebuf);
	}
	/*
	GROUP 1030: Start and End JD, timestep (in JD) per block.
	*/
	fprintf(outfp, spaces); /* blank line */
	sprintf(writebuf, "GROUP   1030");
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]);
	fprintf(outfp, writebuf);
	fprintf(outfp, spaces); /* blank line */

	sprintf(writebuf, "%12.2f%12.2f%11.0f.",
		header->ss[0], header->ss[1], header->ss[2]);
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]); /* pad with spaces */
	fprintf(outfp, writebuf);
	/*
	GROUP 1040: Constant names.
	*/
	fprintf(outfp, spaces); /* blank line */
	sprintf(writebuf, "GROUP   1040");
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]);
	fprintf(outfp, writebuf);
	fprintf(outfp, spaces); /* blank line */

	sprintf(writebuf, "%6d", header->ncon);
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]);
	fprintf(outfp, writebuf);

	/*
	Now write the constant names, 10 per line, each 6 characters long
	preceded by 2 blanks.  Pad names with blanks to make 6 characters.
	*/
	for (i=0; i<header->ncon; i++) {
		fprintf(outfp, "  %-6s", header->cnam[i]);
		if (i % 10 == 9) fprintf(outfp, "\r\n");
	}
	if (i % 10 != 0) {  /* Pad last line with spaces (i is 1 more than above) */
		for ( ; i % 10 != 0; i++) fprintf(outfp, "        ");
		fprintf(outfp, "\r\n");
	}
	/*
	GROUP 1041: Constant values.
	*/
	fprintf(outfp, spaces); /* blank line */
	sprintf(writebuf, "GROUP   1041");
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]);
	fprintf(outfp, writebuf);
	fprintf(outfp, spaces); /* blank line */

	sprintf(writebuf, "%6d", header->nval);
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]);
	fprintf(outfp, writebuf);

	if (header->nval != header->ncon) {
		fprintf(stderr,"Error: number of constants and values not equal.\n\n");
		exit(1);
	}
	/*
	Now read constant values, 3 per line, 26 characters each.
	*/
	for (i = 0; i < header->ncon; i += 3) {
		val1 = header->cval[i];
		val2 = (i+1 < header->ncon) ? header->cval[i+1] : 0.0;
		val3 = (i+2 < header->ncon) ? header->cval[i+2] : 0.0;
		sprintf(writebuf, "%25.17E %25.17E %25.17E   \r\n", val1, val2, val3);
		/* Note that the character holding the sign for each # is left as is. */
		ephcom_doublestrc2f(&writebuf[01]);
		ephcom_doublestrc2f(&writebuf[27]);
		ephcom_doublestrc2f(&writebuf[53]);
		fprintf(outfp, "%s", writebuf);
	}
	/*
	GROUP 1050: Constant values.
	*/
	fprintf(outfp, spaces); /* blank line */
	sprintf(writebuf, "GROUP   1050");
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]);
	fprintf(outfp, writebuf);
	fprintf(outfp, spaces); /* blank line */
	/*
	If there are no coefficients for an ipt[i][] object (i.e., ipt[i][1]==0),
	then ipt[i][0] should contain the value of the next available coefficient
	number rather than 0, as per communication of Myles Standish to Paul Hardy
	on preferred format of ephemeris headers.

	If there are no libration coefficients (i.e., lpt[1]==0), then lpt[0]
	should contain the value of the next available coefficient number rather
	than 0 as well, as per the same communication from Myles Standish.
	*/
	/* First set j to maximum index into ipt[] that has coefficients */
	j = 0;
	for (i=1; i<12; i++)
		if (header->ipt[i][1] > 0 && header->ipt[i][0] > j)
			j = i;
	/* Now set j to next available index count. */
	if (header->lpt[1] > 0 && header->lpt[0] > j)
		j = header->lpt[1] + header->lpt[1] * header->lpt[2] * 3;
	else
		j = header->ipt[j][0] +
		header->ipt[j][1] * header->ipt[j][2] * (j==11 ? 2 : 3);
	for (i=1; i<12; i++)
		if (header->ipt[i][0] == 0) header->ipt[i][0] = j;
	if (header->lpt[0] == 0) header->lpt[0] = j;
	/*
	Write ipt array in transposed order (arrays are transposed in FORTRAN
	from their order in C).
	*/
	for (i = 0; i < 3; i++) {
		for (j = 0; j < 12; j++) {
			fprintf(outfp, "%6d", header->ipt[j][i]);
		}
		fprintf(outfp, "%6d  \r\n", header->lpt[i]);
	}
	/*
	GROUP 1070: Constant values.
	*/
	fprintf(outfp, spaces); /* blank line */
	sprintf(writebuf, "GROUP   1070");
	k = strlen(writebuf);
	strcpy(&writebuf[k], &spaces[k]);
	fprintf(outfp, writebuf);
	fprintf(outfp, spaces); /* blank line */
	/*
	Now we're pointing to the first block of coefficient data, after header.
	*/
	return(0);
}


/*
	Write coefficient block information in ASCII format.
*/
int ephcom_writeascii_block( FILE *outfp, struct ephcom_Header *header, int blocknum, double *datablock) 
{
	double val1, val2, val3; /* To write text line with 3 double precision words */
	int i;

	char writebuf[EPHCOM_MAXLINE + 1];

	/*
	Write first line in block, which is block number and ncoeff.
	Note that lines in the data block files end in "\r\n", while
	lines in the header files just end in "\n".
	*/
	fprintf(outfp, "%6d%6d", blocknum + 1, header->ncoeff);
	for (i=0; i<68; i++) fputc(' ', outfp);
	fprintf(outfp, "\r\n");
	/*
	Now write the data, 3 coefficients per line, 26 characters each.
	Convert format to match what appears in JPL Ephemeris ASCII files.
	*/
	for (i = 0; i < header->ncoeff; i += 3) {
		val1 = datablock[i];
		val2 = (i+1) < header->ncoeff ? datablock[i+1] : 0.0;
		val3 = (i+2) < header->ncoeff ? datablock[i+2] : 0.0;
		/*
		Write values, 3 coefficients per line, pad lines with 0.0E+00
		*/
		sprintf(writebuf,
			"%25.17E %25.17E %25.17E   \r\n",
			val1, val2, val3);
		// printf("[%d]%s==>", strlen(writebuf), writebuf);
		/*
		Now re-format numbers the way the JPL header file writes them:
		all with a leading "0.", so the exponent is one greater.
		Note that here we start at strlen(writebuf)-6, but in the section
		that handles header data coefficients we start at strlen(writebuf)-5.
		This is because the data blocks end ASCII lines with "\r\n" whereas
		the header ends ASCII lines with just "\n".
		*/
		ephcom_doublestrc2f(&writebuf[01]);  /* Reformat first number */
		ephcom_doublestrc2f(&writebuf[27]);  /* Reformat second number */
		ephcom_doublestrc2f(&writebuf[53]);  /* Reformat third number */
		// printf("[%d]%s", strlen(writebuf), writebuf);
		fprintf(outfp, "%s", writebuf);
	}
	return(0);
}

/*
	Print an integer (4-byte) value to the given file with bytes swapped
	if necessary to match network order (Big Endian).  On Intel 80x86
	the bytes will get swapped, on Motorola or SPARC they won't.
*/
int ephcom_outint(FILE *outfp, unsigned u) 
{
	unsigned u2;
	
	fwrite(&u2, 4, 1, outfp);
	return(0);
}

/*
	Print a double precision value to the given file with bytes swapped
	if necessary to match network order (Big Endian).  On Intel 80x86
	the bytes will get swapped, on Motorola or SPARC they won't.
*/
int ephcom_outdouble(FILE *outfp, double x) 
{
	unsigned char ch[8];

	memcpy((void *)ch, (const void *)&x, 8);
	fwrite(ch, 1, 8, outfp);
	return(0);
}

/*
	Write a JPL Ephemeris header in binary format.
*/
int ephcom_writebinary_header(FILE *outfp, struct ephcom_Header *header) 
{
	int blockout; /* number of bytes we've written to current block/rec in file */
	int blockbytes; /* number of bytes in a block, equals 8 * ncoeff */

	int i, j;
	int idate[6];
	char *month[12] = {"JAN","FEB","MAR","APR","MAY","JUN",
		"JUL","AUG","SEP","OCT","NOV","DEC"};

	/*
	Point to beginning of output file.
	*/
	rewind(outfp);
	/*
	First header line: KSIZE= # NCOEFF= #
	*/
	blockbytes = sizeof(double) * header->ncoeff;
	/*
	Start writing output ephemeris, beginning with header.
	*/
	/*
	Header title lines with dates, for example:

	JPL Planetary Ephemeris DE405/LE405
	Start Epoch: JED=  2305424.5 1599 DEC 09 00:00:00
	Final Epoch: JED=  2525008.5 2201 FEB 20 00:00:00
	*/
	sprintf(header->ttl[0],"JPL Planetary Ephemeris DE%03d/LE%03d",
		header->numde, header->numle);
	for (i=strlen(header->ttl[0]); i<84; i++) header->ttl[1][i] = ' ';
	ephcom_jd2cal(header->ss[0], idate, 0);
	sprintf(header->ttl[1],"Start Epoch: JED=%11.1f%5d %3s %02d %02d:%02d:%02d",
		header->ss[0], idate[0], month[idate[1]-1], idate[2],
		idate[3], idate[4], idate[5]);
	for (i=strlen(header->ttl[1]); i<84; i++) header->ttl[1][i] = ' ';
	ephcom_jd2cal(header->ss[1], idate, 0);
	sprintf(header->ttl[2],"Final Epoch: JED=%11.1f%5d %3s %02d %02d:%02d:%02d",
		header->ss[1], idate[0], month[idate[1]-1], idate[2],
		idate[3], idate[4], idate[5]);
	for (i=strlen(header->ttl[2]); i<84; i++) header->ttl[2][i] = ' ';
	header->ttl[0][84] = header->ttl[1][84] = header->ttl[2][84] = '\0';

	/*
	ephcom_Header title lines.

	Write the three title lines to the output file, padded with blanks,
	84 characters long (84 is the first even multiple of 6 that is > 80,
	so the JPL software uses that value because it reads in Fortran 'A6'
	character strings.
	*/
	fprintf(outfp, "%-84s%-84s%-84s", header->ttl[0], header->ttl[1], header->ttl[2]);
	blockout = 3*84;  /* Just wrote 3 84-byte strings to start output file */
	/*
	Now output 400 cnam entries to the output file.
	*/
	for (i = 0; i < header->ncon; i++) {
		fprintf(outfp, "%-6s", header->cnam[i]);
		blockout += 6;
	}
	for ( ; i < 400; i++) {
		fprintf(outfp, "      ");  /* Round out to 400 entries, all blank at end */
		blockout += 6;
	}
	/*
	Binary values: Make sure bytes are in big-endian (network) order for file.
	*/
	for (i=0; i<3; i++) 
	{
		ephcom_outdouble(outfp, header->ss[i]);  /* Write net-order bytes from double precision */
		blockout += 8;
	}
	ephcom_outint(outfp, header->ncon);
	blockout += 4;
	ephcom_outdouble(outfp, header->au);
	blockout += 8;
	ephcom_outdouble(outfp, header->emrat);
	blockout += 8;
	/*
	If there are no coefficients for an ipt[i][] object (i.e., ipt[i][1]==0),
	then ipt[i][0] should contain the value of the next available coefficient
	number rather than 0, as per communication of Myles Standish to Paul Hardy
	on preferred format of ephemeris headers.

	If there are no libration coefficients (i.e., lpt[1]==0), then lpt[0]
	should contain the value of the next available coefficient number rather
	than 0 as well, as per the same communication from Myles Standish.
	*/
	/* First set j to maximum index into ipt[] that has coefficients */
	j = 0;
	for (i=1; i<12; i++)
		if (header->ipt[i][1] > 0 && header->ipt[i][0] > j)
			j = i;
	/* Now set j to next available index count. */
	if (header->lpt[1] > 0 && header->lpt[0] > j)
		j = header->lpt[1] + header->lpt[1] * header->lpt[2] * 3;
	else
		j = header->ipt[j][0] +
		header->ipt[j][1] * header->ipt[j][2] * (j==11 ? 2 : 3);
	for (i=1; i<12; i++)
		if (header->ipt[i][0] == 0) header->ipt[i][0] = j;
	if (header->lpt[0] == 0) header->lpt[0] = j;

	for (j = 0; j < 12; j++) {
		for (i = 0; i < 3; i++) {
			ephcom_outint(outfp, header->ipt[j][i]);
			blockout += 4;
		}
	}
	ephcom_outint(outfp, header->numde);
	blockout += 4;
	for (i = 0; i < 3; i++) {
		ephcom_outint(outfp, header->lpt[i]);
		blockout += 4;
	}
	/*
	Now pad the end of the first record with null bytes.  Note: the
	JPL Fortran software just skips to next record at this point.
	*/
	for (i = blockout; i < blockbytes; i++) {
		fputc('\0', outfp);
	}
	/*
	End of first block.  Now set blockout to 0 and start with next block.
	*/
	blockout = 0;
	for (i=0; i<header->ncon; i++) {
		ephcom_outdouble(outfp, header->cval[i]);
		blockout += 8;
	}
	/*
	Pad with double-precision zeroes for rest of array.
	*/
	for ( ; i < 400; i++) {
		ephcom_outdouble(outfp, (double)0.0);
		blockout += 8;
	}
	/*
	Pad with nulls for rest of block.
	*/
	for (i = blockout; i < blockbytes; i++) {
		fputc('\0', outfp);
	}
	/*
	Finished normally.
	*/
	return(0);
}

/*
	Write a block of data coefficients in JPL binary file format.
*/
int ephcom_writebinary_block( FILE *outfp, struct ephcom_Header *header, int blocknum, double *datablock) 
{
	int i;
	//int ephcom_outdouble(FILE *, double);
	int filebyte;
	int filepos;

	/*
	Find out where we need to point in the binary file.
	*/
	filebyte = (blocknum + 2) * header->ncoeff * 8; /* 8 bytes per coefficient */
	/*
	If the file isn't that large, pad it with null bytes
	*/
	fseek(outfp, 0L, SEEK_END);
	filepos = ftell(outfp);
	if (filepos < filebyte) {
		for (i=0; i < (filebyte - filepos); i++) {
			fputc('\0', outfp);
		}
	}
	/*
	Now go to position where we want to start writing.
	*/
	fseek(outfp, filebyte, SEEK_SET);
	for (i = 0; i < header->ncoeff; i++) {
		ephcom_outdouble(outfp, datablock[i]);
	}

	return(0);
}


/*
	ephcom_parse_block() - Parse a binary block of data.  Warning: verbose!
	Writes parsed output to file pointer outfp.
*/
int ephcom_parse_block(  FILE *outfp, struct ephcom_Header *header, double *datablock) 
{
	int i0, i1, i2, i3;
	int blockword;
	/*
	Names of the objects in Chebyshev coefficient arrays.
	*/
	static char *ephcom_coeffname[13] = {
	"Mercury","Venus","EMBary","Mars","Jupiter","Saturn","Uranus","Neptune",
	"Pluto","Moon","Sun","Nutation","Libration"};

	blockword=0;
	fprintf(outfp, "@%04d StartJD\t%12.2f\n", blockword++, datablock[0]);
	fprintf(outfp, "@%04d EndJD\t%12.2f\n", blockword++, datablock[1]);
	for (i0=0; i0<13; i0++) /* For all bodies */
	{  
		fprintf(outfp, "Body\t%d (%s)\n", i0+1, ephcom_coeffname[i0]);
		for (i1=0; i1<header->ipt[i0][2]; i1++)  /* For all subintervals */
		{
			fprintf(outfp, "  Subinterval %d of %d\n", i1+1, header->ipt[i0][2]);
			for (i2=0; i2 < (i0==11 ? 2 : 3); i2++)   /* For all coordinates */
			{
				fprintf(outfp, "    %cCoefficients\n", 'X' + i2);
				for (i3=0; i3<header->ipt[i0][1]; i3++)   /* For all coefficients */
				{
					blockword =     header->ipt[i0][0] +
						i1*header->ipt[i0][1] * (i0==11 ? 2 : 3) +
						i2*header->ipt[i0][1] + i3 - 1;
					fprintf(outfp, "      @%04d [%2d of %2d] %25.17E\n",
						blockword, i3+1, header->ipt[i0][1], datablock[blockword]);
				}
			}
		}
	}
	return(0);
}


/*
Planetary Ephemeris.  Takes coordinates already calculated in
coords structure an converts to vectors and vector dot in testr[].
Bodies start at 1 for Mercury, to match the JPL PLEPH() numbering.
Values for ntarg and ncntr correspond to locations ntarg-1 and
ncntr-1 in coords->pv[].
*/
int CEphemDE::ephcom_pleph(struct ephcom_Coords *coords, int ntarg, int ncntr, double *r) 
{
	int i;

	if (ntarg != 14 && ncntr != 14) 
	{ /* If not nutation, handle normally */
		if (ntarg == 15 || ncntr == 15) 
		{ /* Libration */
			for (i=0; i<6; i++) 
				r[i] = coords->pv[14][i];
			//    for (i=0; i<6; i++) printf("\nlibration[%d] = %g\n", i, r[i]);
		}
		else 
		{
			for (i=0; i<6; i++) 
			{
				r[i] = coords->pv[ntarg-1][i] - coords->pv[ncntr-1][i];
				//       printf("\n%g = %g - %g\n",
				//              r[i], coords->pv[ntarg-1][i], coords->pv[ncntr-1][i]);
			}
		}
	}
	else 
	{ /* Nutation */
		r[0] = coords->pv[13][0];
		r[1] = coords->pv[13][1];
		r[2] = coords->pv[13][2];
		r[3] = coords->pv[13][3];
		r[4] = 0.0;
		r[5] = 0.0;
		// for (i=0; i<4; i++) printf("\nnutation[%d] = %g\n", i, r[i]);
	}

	return(0);
}

/*
ephcom_cal2jd() - convert calendar date and time to JD.

idate: integer year, month, day, hour, minute, second
calendar_type: -1=Julian; 0=Automatic; 1=Gregorian
return value: double precision Julian Day of idate[]

From pp. 604, 606 in the Explanatory Supplement to the Astronomical Almanac.
*/
double ephcom_cal2jd(int idate[6], int calendar_type) 
{
	double tjd;
	int jd;

	/*
	Convert hours, minutes, seconds to fractional JD.
	*/
	tjd = (idate[3] + (idate[4] + idate[5] / 60.0) / 60.0) / 24.0 - 0.5;
	/*
	Julian calendar.  Explanatory Supplement to Astronomical Alamanac, p. 606.
	If automatic, use Julian calendar for dates before 15 October 1582.
	*/
	if (calendar_type == -1 ||
		(calendar_type == 0 && 
		(idate[0] < 1582 ||                         /* Before 1582 */
		(idate[0] == 1582 &&
		(idate[1] < 10 ||                         /* Before October 1582 */
		(idate[1] == 10 && idate[2] < 15)))))) { /* Before 15 October 1582 */
			jd = 367 * idate[0] -
				(7 * (idate[0] + 5001 + (idate[1] - 9) / 7)) / 4 +
				(275 * idate[1]) / 9 +
				idate[2] + 1729777;
	}
	/*
	Gregorian calendar.
	*/
	else { /* Explanatory Supplement to Astronomical Almanac, p. 604 */
		jd = (1461 * (idate[0] + 4800 + (idate[1] - 14) / 12)) / 4 +
			(367 * (idate[1] - 2 - 12 * ((idate[1] - 14) / 12))) / 12 -
			(3 * ((idate[0] + 4900 + (idate[1] - 14) / 12) / 100)) / 4 +
			idate[2] - 32075;
	}
	/*
	Return value is whole JD number plus fractional JD number.
	*/
	tjd += (double)jd;

	return(tjd);
}

